# R Script analyzes results from conjoint analysis following the Pre Analysis Plan
#
# author: Jonas Markgraf; Guillermo Rosas 
# last update: Mar 23, 2023
###############################################################################

rm(list = ls())
# load packages
library(cjoint)
library(repmis)
library(tidyr)
library(FindIt)
library(readxl)
library(stringr)
library(cregg)
library(sandwich)
library(tidyverse)

# R Version
print (version)
# platform       x86_64-apple-darwin17.0     
# arch           x86_64                      
# os             darwin17.0                  
# system         x86_64, darwin17.0          
# status                                     
# major          4                           
# minor          2.0                         
# year           2022                        
# month          04                          
# day            22                          
# svn rev        82229                       
# language       R                           
# version.string R version 4.2.0 (2022-04-22)
# nickname       Vigorous Calisthenics       
# RStudio        2023.03.0+386 

# set wd
possibles <- c("~/Desktop/replication_files/") # set your wd here

set_valid_wd(possibles)

# load data
surveyData <- readRDS("cleanedConjointData.csv")


## 1. Average Marginal Component Effects --------------------------------------

# specify order of levels
surveyData$consum.crdt <- factor(surveyData$consum.crdt, ordered=TRUE,
  levels=c("5%", "15%","25%"))
surveyData$mrtg.crdt <- factor(surveyData$mrtg.crdt, ordered=TRUE,
  levels=c("2%", "6%","10%"))
surveyData$welfare <- factor(surveyData$welfare, ordered=TRUE,
  levels=c("None", "33% of average income","66% of average income"))
surveyData$unempl <- factor(surveyData$unempl, ordered=TRUE,
  levels=c("0 months", "Up to 12 months","Up to 24 months"))
surveyData$tax <- factor(surveyData$tax, ordered=TRUE,
  levels=c("15% tax", "25% tax","35% tax"))

## Table H.1: --------------------------------------------------
fit.main <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData, 
  pair.id = surveyData$taskId,
  diff=TRUE, 
  cluster=surveyData$ResponseId,
  nway=2)

# refit with test.CausalANOVA
fit.mainTest <- test.CausalANOVA(fit.main,
                                      newdata = surveyData,
                                      diff = T,
                                      pair.id = surveyData$taskId,
                                      cluster=surveyData$id)
summary(fit.mainTest)

# loop to export --------------------------------------

# loop over plots
crdtVars <- c("mrtg.crdt", "consum.crdt")
wlfrVars <- c("welfare", "unempl", "tax")

### THIS IS FIGURE H1 
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/mainANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.mainTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# EFFECT HETEROGENEITY --------------------------------------------------------

# Table H.2: LR scale -----------------------------------

## left
fit.lrLeft <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$lr.cat == "left",], 
  pair.id = surveyData[surveyData$lr.cat == "left",]$taskId,
  cluster=surveyData[surveyData$lr.cat == "left",]$id,
  nway=2)

fit.lrLeftTest <- test.CausalANOVA(fit.lrLeft
  , newdata = surveyData[surveyData$lr.cat == "left",]
  , diff = T, pair.id = surveyData[surveyData$lr.cat == "left",]$taskId
  , cluster=surveyData[surveyData$lr.cat == "left",]$id)

### THIS IS FIGURE H2, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/lrLeftANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.lrLeftTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()


## right
fit.lrRight <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$lr.cat == "right",], 
  pair.id = surveyData[surveyData$lr.cat == "right",]$taskId,
  diff=TRUE, 
  cluster=surveyData[surveyData$lr.cat == "right",]$id,
  nway=2)
fit.lrRightTest <- test.CausalANOVA(fit.lrRight
  , newdata = surveyData[surveyData$lr.cat == "right",]
  , diff = T, pair.id = surveyData[surveyData$lr.cat == "right",]$taskId
  , cluster=surveyData[surveyData$lr.cat == "right",]$id)

### THIS IS FIGURE H2, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/lrRightANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.lrRightTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.3: Welfare Beneficiaries ----------------------

## Beneficiaries
fit.wlfrBenefic <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$wlfrBenefit == 1,], 
  pair.id = surveyData[surveyData$wlfrBenefit == 1,]$taskId,
  diff=TRUE, 
  cluster=surveyData[surveyData$wlfrBenefit == 1,]$id,
  nway=2)
fit.wlfrBeneficTest <- test.CausalANOVA(fit.wlfrBenefic
  , newdata = surveyData[surveyData$wlfrBenefit == 1,]
  , diff = T, pair.id = surveyData[surveyData$wlfrBenefit == 1,]$taskId
  , cluster=surveyData[surveyData$wlfrBenefit == 1,]$id)

### THIS IS FIGURE H3, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/wlfrBeneficANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.wlfrBeneficTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()


## No Beneficiaries
fit.wlfrNoBenefic <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$wlfrBenefit == 0,], 
  pair.id = surveyData[surveyData$wlfrBenefit == 0,]$taskId,
  diff=TRUE, 
  cluster=surveyData[surveyData$wlfrBenefit == 0,]$id,
  nway=2)
fit.wlfrNoBeneficTest <- test.CausalANOVA(fit.wlfrNoBenefic
  , newdata = surveyData[surveyData$wlfrBenefit == 0,]
  , diff = T, pair.id = surveyData[surveyData$wlfrBenefit == 0,]$taskId
  , cluster=surveyData[surveyData$wlfrBenefit == 0,]$id)

### THIS IS FIGURE H3, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/wlfrNoBeneficANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.wlfrNoBeneficTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.4: Precarious Workers -------------------------

## No Precarious Employment
fit.noPrecEmpl <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$precEmpl == 0,], 
  pair.id = surveyData[surveyData$precEmpl == 0,]$taskId,
  diff=TRUE, 
  cluster = surveyData[surveyData$precEmpl == 0,]$id,
  nway=2)

fit.noPrecEmplTest <- test.CausalANOVA(fit.noPrecEmpl
  , newdata = surveyData[surveyData$precEmpl == 0,]
  , diff = T, 
  pair.id = surveyData[surveyData$precEmpl == 0,]$taskId
  , cluster=surveyData[surveyData$Q61 == 0,]$id)

### THIS IS FIGURE H4, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/noPrecEmplANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.noPrecEmplTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Precarious Employment
fit.PrecEmpl <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$precEmpl == 1,], 
  pair.id = surveyData[surveyData$precEmpl == 1,]$taskId,
  diff=TRUE, 
cluster= surveyData[surveyData$precEmpl == 1,]$id,
  nway=2)

fit.PrecEmplTest <- test.CausalANOVA(fit.PrecEmpl
  , newdata = surveyData[surveyData$precEmpl == 1,]
  , diff = T, pair.id = surveyData[surveyData$precEmpl == 1,]$taskId
  , cluster=surveyData[surveyData$precEmpl == 1,]$id)

### THIS IS FIGURE H4, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/PrecEmplANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.PrecEmplTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.5: Unemployment Risk --------------------------

# High Unemployment Risk 
fit.HiRisk <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$unempRisk3.cat == "High Risk",], 
  pair.id = surveyData[surveyData$unempRisk3.cat == "High Risk",]$taskId,
  diff=TRUE, 
  cluster= surveyData[surveyData$unempRisk3.cat == "High Risk",]$id,
  nway=2)
fit.HiRiskTest <- test.CausalANOVA(fit.HiRisk
  , newdata = surveyData[surveyData$unempRisk3.cat == "High Risk",]
  , diff = T, pair.id = surveyData[surveyData$unempRisk3.cat == "High Risk",]$taskId
  , cluster=surveyData[surveyData$unempRisk3.cat == "High Risk",]$id)

### THIS IS FIGURE H5, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/riskHiANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.HiRiskTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Low Unemployment Risk
fit.LoRisk <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$unempRisk3.cat == "Low Risk",], 
  pair.id = surveyData[surveyData$unempRisk3.cat == "Low Risk",]$taskId,
  diff=TRUE, 
  cluster= surveyData[surveyData$unempRisk3.cat == "Low Risk",]$id,
  nway=2)
fit.LoRiskTest <- test.CausalANOVA(fit.LoRisk
  , newdata = surveyData[surveyData$unempRisk3.cat == "Low Risk",]
  , diff = T, pair.id = surveyData[surveyData$unempRisk3.cat == "Low Risk",]$taskId
  , cluster=surveyData[surveyData$unempRisk3.cat == "Low Risk",]$id)

### THIS IS FIGURE H5, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/riskLoANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.LoRiskTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.6: Sectoral Unemployment ---------------------------
 
# High Unemployment Risk 
fit.HiSUR <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$highUnemp == "High",], 
  pair.id = surveyData[surveyData$highUnemp == "High",]$taskId,
  diff=TRUE, 
  cluster= surveyData[surveyData$highUnemp == "High",]$id,
  nway=2)
fit.HiSURTest <- test.CausalANOVA(fit.HiSUR
  , newdata = surveyData[surveyData$highUnemp == "High",]
  , diff = T, pair.id = surveyData[surveyData$highUnemp == "High",]$taskId
  , cluster=surveyData[surveyData$highUnemp == "High",]$id)

### THIS IS FIGURE H6, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/SURHiANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.HiSURTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Low Unemployment Risk
fit.LoSUR <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data=surveyData[surveyData$highUnemp == "Low",], 
  pair.id = surveyData[surveyData$highUnemp == "Low",]$taskId,
  diff=TRUE, 
  cluster= surveyData[surveyData$highUnemp == "Low",]$id,
  nway=2)
fit.LoSURTest <- test.CausalANOVA(fit.LoSUR
  , newdata = surveyData[surveyData$highUnemp == "Low",]
  , diff = T, pair.id = surveyData[surveyData$highUnemp == "Low",]$taskId
  , cluster=surveyData[surveyData$highUnemp == "Low",]$id)

### THIS IS FIGURE H6, LEFT 
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/SURLoANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.LoSURTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.7: Difficulty to obtain loan ------------------
surveyData$crdtAccess2 <- ifelse(surveyData$crdtAccess >= 4 & !is.na(surveyData$crdtAccess)
  , "Difficult", ifelse(surveyData$crdtAccess <= 2 & !is.na(surveyData$crdtAccess)
    , "Easy", "Else"))

# Easy Credit Access
fit.crdtAcsEasy <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$crdtAccess2 == "Easy",], 
  pair.id = surveyData[surveyData$crdtAccess2 == "Easy",]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$crdtAccess2 == "Easy",]$id,
  nway=2)
fit.crdtAcsEasyTest <- test.CausalANOVA(fit.crdtAcsEasy
  , newdata = surveyData[surveyData$crdtAccess2 == "Easy",]
  , diff = T, pair.id = surveyData[surveyData$crdtAccess2 == "Easy",]$taskId
  , cluster=surveyData[surveyData$crdtAccess2 == "Easy",]$id)

### THIS IS FIGURE H7, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/crdtAcsEasyANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.crdtAcsEasyTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Difficult Credit Access
fit.crdtAcsDifficult <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$crdtAccess2 == "Difficult",], 
  pair.id = surveyData[surveyData$crdtAccess2 == "Difficult",]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$crdtAccess2 == "Difficult",]$id,
  nway=2)
fit.crdtAcsDifficultTest <- test.CausalANOVA(fit.crdtAcsDifficult
  , newdata = surveyData[surveyData$crdtAccess2 == "Difficult",]
  , diff = T, pair.id = surveyData[surveyData$crdtAccess2 == "Difficult",]$taskId
  , cluster=surveyData[surveyData$crdtAccess2 == "Difficult",]$id)

### THIS IS FIGURE H7, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/crdtAcsDifficultANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.crdtAcsDifficultTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.8: Financial Crisis ---------------------------

# Experienced Crisis Losses
fit.crisisLoss <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$expFinCrisis.cat == "Losses in financial crisis",], 
  pair.id = surveyData[surveyData$expFinCrisis.cat == "Losses in financial crisis",]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$expFinCrisis.cat == "Losses in financial crisis",]$id,
  nway=2)
fit.crisisLossTest <- test.CausalANOVA(fit.crisisLoss
  , newdata = surveyData[surveyData$expFinCrisis.cat == "Losses in financial crisis",]
  , diff = T, pair.id = surveyData[surveyData$expFinCrisis.cat == "Losses in financial crisis",]$taskId
  , cluster=surveyData[surveyData$expFinCrisis.cat == "Losses in financial crisis",]$id)

### THIS IS FIGURE H8, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/crisisLossANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.crisisLossTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Experienced no Crisis Losses
fit.crisisNoLoss <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$expFinCrisis.cat == "No losses in financial crisis",], 
  pair.id = surveyData[surveyData$expFinCrisis.cat == "No losses in financial crisis",]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$expFinCrisis.cat == "No losses in financial crisis",]$id,
  nway=2)
fit.crisisNoLossTest <- test.CausalANOVA(fit.crisisNoLoss
  , newdata = surveyData[surveyData$expFinCrisis.cat == "No losses in financial crisis",]
  , diff = T, pair.id = surveyData[surveyData$expFinCrisis.cat == "No losses in financial crisis",]$taskId
  , cluster=surveyData[surveyData$expFinCrisis.cat == "No losses in financial crisis",]$id)

### THIS IS FIGURE H8, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/crisisNoLossANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.crisisNoLossTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.9: Income Effect ------------------------------

# High Income
fit.incomeHi <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$income.cat == "High Income",], 
  pair.id = surveyData[surveyData$income.cat == "High Income",]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$expFinCrisis.cat == "High Income",]$id,
  nway=2)
fit.incomeHiTest <- test.CausalANOVA(fit.incomeHi
  , newdata = surveyData[surveyData$income.cat == "High Income",]
  , diff = T, pair.id = surveyData[surveyData$income.cat == "High Income",]$taskId
  , cluster=surveyData[surveyData$income.cat == "High Income",]$id)

### THIS IS FIGURE H9, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/incomeHiANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.incomeHiTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Low Income
fit.incomeLo <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$income.cat == "Low Income",], 
  pair.id = surveyData[surveyData$income.cat == "Low Income",]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$expFinCrisis.cat == "Low Income",]$id,
  nway=2)
fit.incomeLoTest <- test.CausalANOVA(fit.incomeLo
  , newdata = surveyData[surveyData$income.cat == "Low Income",]
  , diff = T, pair.id = surveyData[surveyData$income.cat == "Low Income",]$taskId
  , cluster=surveyData[surveyData$income.cat == "Low Income",]$id)

### THIS IS FIGURE H9, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/incomeLoANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.incomeLoTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.10: Education Effect ---------------------------

# High Education
fit.educHi <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$educ.cat == "Higher Education",], 
  pair.id = surveyData[surveyData$educ.cat == "Higher Education",]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$educ.cat == "Higher Education",]$id,
  nway=2)
fit.educHiTest <- test.CausalANOVA(fit.educHi
  , newdata = surveyData[surveyData$educ.cat == "Higher Education",]
  , diff = T, pair.id = surveyData[surveyData$educ.cat == "Higher Education",]$taskId
  , cluster=surveyData[surveyData$educ.cat == "Higher Education",]$id)

### THIS IS FIGURE H10, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/educHiANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.educHiTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# High Education
fit.educLo <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$educ.cat == "No Higher Education",], 
  pair.id = surveyData[surveyData$educ.cat == "No Higher Education",]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$educ.cat == "No Higher Education",]$id,
  nway=2)
fit.educLoTest <- test.CausalANOVA(fit.educLo
  , newdata = surveyData[surveyData$educ.cat == "No Higher Education",]
  , diff = T, pair.id = surveyData[surveyData$educ.cat == "No Higher Education",]$taskId
  , cluster=surveyData[surveyData$educ.cat == "No Higher Education",]$id)

### THIS IS FIGURE H10, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/educLoANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.educLoTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Table H.11: Homeownership ------------------------------

# Home owner
fit.homeOwner <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$howner == 1,], 
  pair.id = surveyData[surveyData$howner == 1,]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$howner == 1,]$id,
  nway=2)
fit.homeOwnerTest <- test.CausalANOVA(fit.homeOwner
  , newdata = surveyData[surveyData$howner == 1,]
  , diff = T, pair.id = surveyData[surveyData$howner == 1,]$taskId
  , cluster=surveyData[surveyData$howner == 1,]$id)

### THIS IS FIGURE H11, LEFT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/homeOwnerANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.homeOwnerTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()

# Home renter
fit.homeRenter <- CausalANOVA(formula=forcedChoice ~ consum.crdt
  + mrtg.crdt + welfare + unempl + tax,
  data = surveyData[surveyData$howner == 0,], 
  pair.id = surveyData[surveyData$howner == 0,]$taskId,
  diff=TRUE,
  cluster= surveyData[surveyData$howner == 0,]$id,
  nway=2)
fit.homeRenterTest <- test.CausalANOVA(fit.homeRenter
  , newdata = surveyData[surveyData$howner == 0,]
  , diff = T, pair.id = surveyData[surveyData$howner == 0,]$taskId
  , cluster=surveyData[surveyData$howner == 0,]$id)

### THIS IS FIGURE H11, RIGHT
# pdf(file='~/Dropbox/conjointCredit/03_JOPsubmission/replication_files/graphPlots/homeRenterANOVA_plot.pdf', width=10, height=12)
par(mfrow = c(3, 2), mar = c(1,1,1,1))
for(j in wlfrVars) {
  for(i in crdtVars) {
    print(plot(fit.homeRenterTest, type="ConditionalEffect", fac.name=c(j,i))) 
  }
}
# dev.off()
